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Abstract 

We consider the localization problem of multiple wideband sources in a multi-path environment by coherently 
taking into account the attenuation characteristics and the time delays in the reception of the signal. Our 
proposed method leaves the space for unavailability of an accurate signal attenuation model in the environment by 
considering the model as an unknown function with reasonable prior assumptions about its functional space. Such 
approach is capable of enhancing the localization performance compared to only utilizing the signal attenuation 
information or the time delays. In this paper, the localization problem is modeled as a cost function in terms 
of the source locations, attenuation model parameters and the multi-path parameters. To globally perform the 
minimization, we propose a hybrid algorithm combining the differential evolution algorithm with the Levenberg- 
Marquardt algorithm. Besides the proposed combination of optimization schemes, supporting the technical details 
such as closed forms of cost function sensitivity matrices are provided. Finally, the validity of the proposed method 
is examined in several localization scenarios, taking into account the noise in the environment, the multi-path 
phenomenon and considering the sensors not being synchronized. 



1 Introduction 

A challenging and highly demanding signal process- 
ing application is the localization of signal sources 
using the physical measurements at some sensors in 
the environment. Source localization has become 
an important task in various applications such as 
mobile communications, global positioning system 
(GPS), radar, sonar, navigation, seismology and geo- 
physics [l}{5]- 

During the recent decades various algorithms 
have been proposed to estimate the location of the 
signal sources. These methods utilize different sig- 
nal characteristics at different sensors and generally 
can be classified in three main categories: using the 



time difference of arrival (TDOA) ; analyzing the sig- 
nal direction of arrival (DOA) at distinct arrays; 
and using the differences in the signal amplitude or 
received energy level. For a constant propagation 
speed, the TDOA among different sensors is propor- 
tional to the source-sensor range differences and may 
be estimated through methods such as cross corre- 
lation (CC) [6] or its generalized version (GCC) [7]. 
The source locations can then be estimated using 
geometric methods such as linear, spherical or hy- 
perbolic intersections [8 10 . To estimate the DOA, 
for narrowband signals, high resolution algorithms 
such as multiple signal classification (MUSIC) [II] 
and maximum likelihood (ML) |l2] are proposed. 
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In 13 , the authors propose an approximate maxi- 



mum likelihood method (AML) for wideband signals 
using spectral properties of the signal when rather 
long sample streams are available. In this method, 
the corresponding cost function can be directly ex- 
pressed in terms of the source locations or in a far 
field case may be expressed in terms of the relative 
time delays followed by a post processing step to find 
the source locations from the corresponding DOAs. 
The post processing step may be carried out through 
geometric methods such as cross bearing or a ma- 
chine learning approach such as the support vector 
machine (SVM) method 14 . Using the differences 



in the signal intensity or energy level for the purpose 
of localization is a more recent technique 15 16 



Theoretically, this class of localization can be con- 
sidered for both narrowband and wideband signals 
by only taking into account the attenuation infor- 
mation and usually neglecting the time delay infor- 
mation. For these methods, a precise attenuation 
model in the environment is inevitable for an ac- 
curate localization. Moreover, from an optimization 
perspective the resulting cost functions in these kind 
of approaches usually undergo many local optima 
and saddle points which require considering specific 



optimization schemes 17 



In this paper, we tackle the problem of local- 
ization of multiple wideband sources by coherently 
taking into account the TDOA and the amplitude at- 
tenuation pattern. Our method generalizes the AML 
approach to utilize the signal attenuation character- 
istics. We provide a more robust algorithm in which 
the targets should simultaneously satisfy the correct 
time delays among the sensors and provide sensi- 
ble level of attenuation at each sensor. Unlike the 
aforementioned energy and intensity based methods 
which not only ignore the time delay stamps but 
also require knowing the signal attenuation model, 
we benefit using the delay information and as a gen- 
eralization to our recent work in [18] , leave the space 
for not knowing an exact signal attenuation model in 
the environment by suggesting an appropriate func- 
tional space for it. We minimize a cost function 
which is obtained through maximum likelihood ap- 
proach from which the locations, attenuation model 
parameters and the multi-path parameters are ob- 
tained. To apply the minimization we propose a hy- 
brid approach combining the differential evolution 



algorithm 19 with the Levenberg-Marquardt algo- 
rithm 20 . This combination provides a minimiza- 



tion scheme which is likely to globally search for the 



optima and rather quickly converges to the accu- 
rate results. Through simulations and Cramer-Rao 
bound we verify the effectiveness of the novel method 
introduced in this paper. 

This paper is organized as follows. In section 
[2] we propose a general form for the received sig- 
nal at every sensor and later provide an adaptive 
model for the signal attenuation based on Laurent 
polynomials. In section [3j a maximum likelihood es- 
timation of the source location and attenuation pa- 
rameters is proposed. We also provide the Cramer- 
Rao bound for this estimation problem. For the pur- 
pose of minimization in section [4] a hybrid approach 
combining the differential evolution algorithm with 
the Levenberg-Marquardt is proposed for which the 
combination algorithm and closed form equations for 
calculation of the Jacobian are provided. In section 
[5] we examine the efficiency of proposed method 
through some examples and finally there arc some 
concluding remarks in section [6] 



2 Problem Definition 
2.1 Signal Model 

Although the general approach proposed in this pa- 
per is in theory independent of signal nature and the 
type of sensors used, in order to make reasonable 
simulations we consider acoustic source localization. 
Consider N acoustic sources having unknown loca- 
tions rs n - Each source is omni-directionally emit- 
ting a signal s„(i), n = 1,---,N at the time frame t. 
We also consider M acoustic microphones that are 
placed in known positions rM m i m = 1j"">-^j m the 
same environment. For every source in the environ- 
ment, the function that describes the signal atten- 
uation at a specific point is a(p), where p is the 
distance from the point to the source. In general, 
the signal attenuation function may be a function of 
various parameters such as signal frequency, medium 
inhomogeneity, etc. To simplify the problem, in this 
paper we consider this function to be an identical 
form for all sources and solely function of the dis- 
tance to the source. However, unlike some previous 
energy based localizations (e.g., see 16p"7 ) in which 
the attenuation is known to be proportional to 
the actual form of a(-) is considered unknown func- 
tion here. This type of modeling provides an ad- 
ditional flexibility to the problem for more realistic 
scenarios where the inverse proportionality of a(-) to 
p is violated due to other parameters, such as signal 
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bandwidth and medium inhomogeneity. Considering 
s n (t - 1 x N s /v) to be the signal measured 1 length 
unit away from every source with N s being the sam- 
ples per second and v being the propagation speed, 
ideally the overall received signal samples from the 
acoustic sources at every microphone is modeled as 



N 



i)s n (t- 



(1) 



for 



t = 0,l,—,nt-l, m = l,2,—,M. 



',n ll ~ \\ r M„ 

,th 



rs n || is the distance from n 



th 



Here p. 

source to microphone and r TOi „ = p m ^ n N s /v is 
the corresponding time samples delay in receiving 
the signal. The received signal in ([lj is normalized 
to each microphone gain level in order to decrease 
the number of notations. A more realistic model 
takes into account the noise in the environment and 
also the signals going through a multi-path channel 
before arriving at every sensor, hence we rewrite the 
received signal as 

N 

^-m(^) — ^(Pm,n ) — ^*m,n) 

n=\ 

N P m ,n 
71=1 p=l 



+ W m {t). 



(2) 



The term w m (t) represents the background noise 
which is considered to be a zero-mean white Gaus- 
sian with variance a 2 for the sake of this paper; 
i.e., Gaussianity is not a limiting assumption in this 
paper. Between the rfi 1 source and rrfi 1 micro- 
phone we consider P m>ra indirect paths each causing 
7m,n,p loss in the signal amplitude and T m n p de- 
lay in the signal reception, modeling the multi-path 
phenomenon. Beside the positions rs n , which are 
the main unknowns of the localization problem, the 
signals s n (t), the multi-path parameters "f m ,n,p and 
T m .n,p, and the propagation loss function a(-) are 
also unknown and should be determined based on 
the received signals at the sensors. The appearance 
of T m .„ (which is related to the unknown quanti- 
ties rs n ) and T mnp as the argument of an unknown 
signal s n (t) causes an extra complexity for any opti- 
mization scheme performed to solve the localization 
problem. However this problem may be remedied 
by applying the discrete Fourier transform to ^ to 
extract the delays and form an equivalent equation 



in which the unknown parameters are separated in 
individual terms, i.e., 

N j'2tt 
X mU) = T, a(Pm,«)exp( fr m , n )S n (f) 

t7=1 n f 

+ 2^ ^ 7m,n,p ex p( /^m,n,p)^n(/) 
n=l p=l n f 

+ CM) (3) 



for 



/ = 0,l,-,n / -l, 



l,2,-,M. 



Here, X m (f), S n (f) and £ m (f) are the data, signal 



and noise spectrums respectively. As stated in 13 
we emphasize on the fact that for (|3| to be a valid 
equivalent form of ^ , we need n t to be large enough 
to avoid edge effects and accordingly nf > n t . In 
general having more samples from the signal better 
poses the problem. 



2.2 A Low Order Representation of Signal At- 
tenuation Model 

As discussed earlier, our assumption about the at- 
tenuation model in the environment in this paper is 
an identical model for all sources, which is only de- 
pendent on the distance of the point to the acoustic 
source. In an ideal environment, a(p) can be well 
modeled as a multiple of p~ x . Since there are dif- 
ferent parameters involved in the attenuation mod- 
eling, a(p) is being considered as an unknown here. 
However, in order to keep the well-posedness of the 
problem, we choose it to be an element of a low di- 
mensional function space. For this sake, we consider 
a(p) to be a Laurent polynomial of limited order as 



L>0. 



(4) 



In this model, only negative powers of p are consid- 
ered, which is due to the fact that for an attenuation 
model we are physically required to have 



lim a(p) = 0. 



(5) 



In many applications the low order representation 
of a{p) in Q is acceptable enough to model the at- 
tenuation and usually considering only few terms in 
the series (i.e., L rather small), would suffice for the 
localization problem. 
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3 A Maximum Likelihood Estimation of 

the Unknowns 
3.1 Derivation 

Based on the general attenuation model proposed, 
matching of the data spectrum with the model can 
be expressed by using @ in as 

N L '27T 

x m(f) = y y /3<p m ,„exp( fr m ,n)S n (f) 

— n f 



n=l 1=1 

N Pw. 



' Z Z 7™ : n, p exp( ffm,n,p)S n (f) 

71=1 p=l U f 



(6) 



The central limit theorem states that £, m (f), 
which is a transformed zero mean Gaussian ran- 
dom variable to the frequency domain, should 
be a complex zero mean Gaussian with variance 
n t (T 2 . For every frequency bin / having X(f) = 
[X 1 (f),-,X M (f)] T , S(f) = [Si(f),-,S N (f)] T 
and £(/) = [&(/), -,£m(/)] T , (§ can be written 
in a matrix form as 

X(f) = (K(f) + H(f))S(f) + t(f) (7) 

where K(f) = R(f)(3 with 

/3= [Pu-,Pl] t 9 In*n, (8) 

for which <g> represents the Kronecker product and 
InxN the identity matrix of size N x N, and 



H(/) = [fli(/),-,fli(/)], 



(9) 



where 



Re(f) 



Pl \e n r 



U>M,i e " /l 



L /pi 



l /pm,i 



PM,N e " /l 



/Pi 



f Pm,n 



for £ = 1,---,L. The matrix H(f) is related to the 
multi-path parameters and its elements are 

[ if (/)](m,n) = E 7m,n, P exp( /f m ,„, p ). (10) 

Rewriting the negative log-likelihood function to 
estimate the unknown parameters 9 including the 
source positions, source signal spcctrums, multi-path 
parameters and quantities j3g, we have 



9* = argminQ' H Q 



(11) 



where 



Q 



Q(o) 

Q(n f /2) 



(12) 



and Q(f) = X(f) - K(f)S(f) using the notation 

K(f) = K(f)+H(f). (13) 

Similar to the idea in [13], for a real valued signal, 
we can only consider up to n//2 frequency bins and 
form Q with blocks of Q(f) for / = 0, l,-,n//2. We 
would like to highlight the fact that in [13] , the zero 
frequency bin is ignored due to producing a constant 
term in the likelihood function, however in our ap- 
proach the matrices K(0) and H(Q) are still depen- 
dent on p m ^ n and the multi-path parameters 7 m! n.p 
and hence worth being considered. 



Clearly, the minimization in (11) is equivalent 
to minimizing Q H (f)Q(f) for every /. Consider- 
ing the unknown signal spectrum S(f), the minima 
should satisfy 



d{Q H U)Q{f)) 
dS H {f) 



(14) 



which results in S(f) = k\f)X{f) with k\f) 
representing the pseudo-inverse of K(f). Replacing 
the obtained S(f) in Q(f) results in 



Q(f) = x(f)-k(f)k\f)x(f), 



(15) 



for / = 0, Tif/2, and therefore the unknowns of the 
minimization reduce to the source positions, multi- 
path parameters and the attenuation coefficients. 
Considering a 2D localization problem, as the case 
in the example section, neglecting the multi-path the 
vector of unknowns would be 



= [xs 1 ,—,xs H ,ys 1 ,—,Vs K ,0u—,0L]' i 



(16) 



where x$ n and y$ n are the x and y components of 
the position vector r$ n ■ In case of multi-path, the 



parameters 7 



and 



are also included in 0. 



The approach is clearly not only limited to 2D Carte- 
sian systems and 3D Scenarios and other coordinate 
systems may be considered. 
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3.2 Cramer-Rao Lower Bounds for the Estimated 
Parameters 

For an unbiased parameter estimation problem, the 
Cramer-Rao Lower Bound (CLRB) is a theoretical 
lower bound on the variance of the problem esti- 
mates. Based on (JT]), the total model relating the 
parameters of interest to the complete data set is 



X = g(0;S)+S. 



(17) 



Here X = [X(0) T ,-,X(n f /2) T ] T is the full 
data set, S = [S(0) T ,-, S(n f /2) T ] T contains 
the signal spectrum of all the sources and £ = 
[£(0) T ,--,£(7J//2) T ] T is the corresponding noise 
vector. Moreover, Q{9\ S) = KS for which the ma- 
trix K explicitly dependent on 9 is 



K 



K(0) 

o k{\) 














K(n f /2) J 



(18) 



The CRLB is defined as the diagonal elements of the 



inverse Fisher matrix F, which for the model in ( 17) 
is representable as [21] 



Here 



19 : 



S 
9 



(19) 
(20) 



and 7?.^ is the noise covariance matrix which for our 
problem is simply ritO~ 2 I. The matrix [<9<?/9i9] is 
composed of the sub-blocks [dQ/dS], [dQ/drs n ], 
[dQldf3f\, [dg/d-f m , n , p ] and [dg/df m ^ iP ]. Clearly 



dg 
as 



K. 



(21) 



For the 9 parameters, since K is composed of K(f), 
we only discuss the sensitivity of K(f) to every 
class of parameters. Based on the fact that K(f) = 
R(f)(3 we can write 



9KU) dKjf) R(f) dP 

df3 e dp t U) df3 e 



1,2,-,L (22) 



where 



d/3 



£ clement 



= [0,-,0, 1,0,-,0] T ®7 



NxN- 



(23) 



To calculate dK(f)/dxs n , we have 

dk{f) dR(f) 



dxs dx 



(3, n=l,2,-,N (24) 



where 



dR(f) _r dR 1 (f) _ m L (fh (25) 

dx c 



dxs 



dxt 



The matrix dRi(f)ldxs n is a matrix the same size 

as Ri(f), with all columns being zero except the rfi 1 
column. Simply applying the derivative shows that 
the (m,n) element of dR((f)/dxs n is related to the 
(m,n) element of Ri(f) through 



where 



— = <5(n,n) 

n r J rn,n \ pm,n } 



pm,n 

£ j2nN. 



5(n, n') 



Pm.n flf 

1, n = n' 



—) [R*(f)] , 



(26) 



0, n + n' 

An analogous technique is used to derive 

dk{f)id VSn . 

For the multi-path parameters we have 
dk(f)/d~f m . n , p = dH(f)/djm tn . p and also have 
dK(f)/dT m , n , p = dH(f)/dT m ^ hp . Accordingly the 
elements of each matrix are obtained through 



T dH(f) -i , , jztt 
- = b{m,m )d{n,n ) exp( JT m ,n, P ) 

(27) 

and 

T dH(f) ] , , j2n 

-— = -d{m,m )d{n,n ) /7m,n,p 

l OT m < „'„J(m,n) n f 



*e*p(-^f„,„,). 

n f 



(28) 



Specifying the elements of the Fisher matrix F yields 
the CRLB values for all the estimations. 



4 Minimization Strategy 

The minimization in ( |11[ ) may be performed through 
various optimization schemes, most generally cat- 
egorized as global and local optimizations. For 
a global optimization different approaches such 
as deterministic, stochastic or evolutionary and 
metaheuristic methods may be considered [22}|24] . 
Clearly for an accurate localization, global minimiz- 



ers of (11) are required. However in general, us- 



ing global methods to optimize an arbitrary function 



5 



may be iteratively or computationally expensive. As 
an alternative to this and specifically for a least 



squares cost function as (11), local search methods 



such as gradient descent and quasi-Newton methods 
may be considered 
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Although these methods 
can be relatively faster than the global ones, there 
is always a chance of getting trapped into a local 
minima. In the context of localization, although for 
good initial estimates of the source relatively fast 
methods such as the gradient descent and alternat- 
ing projection are proposed, to increase the chances 
of finding a global minima the process usually in- 
volves exhaustive search methods such as the grid 



13 16 



search and multiresolution search 

For the purpose of this paper we consider a hy- 
brid approach combining a global search method 
with a fast local search method 25 26 . Hybrid 



methods have received considerable interests in dif- 

More specif- 



ferent areas in the recent years [26 - 29 



ically we consider a hybrid combination of the Dif- 
ferential Evolution algorithm (DE) [l9j as successful 
evolutionary search with the Levenberg-Marquardt 



algorithm (LMA) 20 30 as a rather fast and robust 



local search method. Before getting to the combina- 
tion scheme, we provide a brief description of each 
method highlighting the main technical issues specif- 
ically in the context of our localization problem. 



4.1 Differential Evolution Algorithm 

DE is among the metaheuristic and evolutionary 
global optimization schemes. Simplicity and suc- 
cessful performance are the main advantages of this 
algorithm. Considering 6 = [9i,02,-",0d] to be the 
vector of problem unknowns with size D, at every 
generation G of the algorithm Np parameter vec- 
tors 9 i:G = [0i,i,G,02,i,G,---,OD,i,G], i = 1,2,— ,N P , 
are generated. The initial population is randomly 
chosen with a uniform distribution in the search re- 
gion. For this work we consider the DE/rand/l/bin, 
which is a general and widely used strategy of this 
For every generation three main 



algorithm 19 31 



operations are performed as follows. 
4-1- 1 Mutation 

In this phase a mutant vector /x, G is generated as 

Mi)G = ri , G + F(0 r2>G -0 r3 , G ), (29) 

where r\, T2 and r% are randomly selected indices 
among 1,2, —,Np and F e [0,2] is a constant real 



scalar controlling the difference vector amplification. 



4-1-2 Crossover 

A mixing with the mutant vector to increase the di- 
versity of the population is performed by generating 
new trial vectors v iG of length D, defined as 



VdA,G 



u d,i,G, 
Qd,i,G, 



r(d) [0! i] < C R ord= k(i) 
otherwise, 



(30) 

with d = 1,2, D. Here Cr e [0, 1] is the crossover 
constant, r(<i)[o.i] is the d^ 1 evaluation of a uni- 
form random number generator in [0, 1] and k(i) e 
{1, 2, D} is a randomly chosen index ensuring that 
Vi : G takes at least one of the elements of \i i G . 



4-1-3 Selection 

At this step a next generation population member 
Oi^G+i is produced by a selection among 9iG and 
Vi : G- This selection is based on the fitness, and ba- 
sically, the vector with the lower cost proceeds to the 
next generation. 



4.2 A Levenberg-Marquardt Algorithm for the 
local Minimization 

As the local minimization scheme, we suggest using 
the LMA. Our attention towards this algorithm is 
based on several advantages. LMA is basically con- 
sidered as a Newton type method and provides a 
rather quadratic convergence. Meanwhile this algo- 
rithm benefits from stability and uses a trust region 
approach 



30 



The other feature of this method, 
considered as an advantage over other methods such 
as the gradient descent, is its suitability for cases 
where there are different variables of different types 
as the cost function arguments. In fact LMA is al- 
most independent of variable scaling, while for meth- 
ods such as the gradient descent, minimizing a cost 
function dependent on a set of variables with differ- 
ent natures and scales requires appropriate param- 



eter scaling to guarantee a proper convergence 30 
This is a demanding feature for our problem where 
the 8 vector in general consists of the source loca- 
tions, attenuations coefficients and the multi-path 
parameters. 

In the LMA which is an iterative algorithm, we 
start with a 9^ as the starting point. At every it- 
eration, having 9^' already in hand, can be 
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obtained by solving 

{Je T Je + A^I)^" 1 ) - 0«) = -J S T Q, (31) 

where Q is the vertical vector of length Mrif/2 
shown in ( 12 ) and obtained for values 9^' at that 
iteration. The parameter X^' is the damping fac- 
tor, obtained at every iteration based on the trust 
region approach 20 30 . The Jacobian matrix Jg 



contains the sensitivities of Q to every element of 0. 
In order to run the algorithm we need to know the 
Jacobian matrix at every iteration, obtaining which 
is discussed in the Appendix. 



4.3 The Hybrid Combination Scheme 

For the purpose of combining the DE with the LMA, 
we propose using a sequential hybridization approach 



26 



In this approach, the DE initially starts the 
minimization by generating consecutively more fit- 
ting generations. After a certain number of gener- 
ations or after getting relatively slow in decreasing 
the fitness, the best in the last generation is passed 
to the LMA algorithm as an initialization. The min- 
imization continues until convergence. An illustra- 
tion of this algorithm is provided in Fig. [l] 



5 Simulation Results 

To examine the method developed in the previ- 
ous section, we consider some localization examples 
in this section. In the first example we consider 
a reverberation-free environment to show the effi- 
ciency of the method for such cases and provide a 
comparative study for this scenario. The second ex- 
ample brings more realistic issues such as the multi- 
path, and sensor synchronization error into the prob- 
lem and examines the performance of the proposed 
method for such cases. 

Before proceeding with the examples we would 
like to highlight a fact regarding the relationship be- 
tween the cost function and the matrix K . Referring 
to (15), it can be easily verified that scaling K(f) 
by a scalar does not change the cost function. In 



other words if the (3i and 7 m ,n,p values are simulta- 
neously scaled by a scalar value, the cost function 
remains the same. Therefore, we rewrite the atten- 
uation model in Q as 



which somehow normalizes a(-) and unifies the rep- 
resentation. Clearly, since the desired unknowns of 
the problem are the acoustic source coordinates, ob- 
taining a multiple of the attenuation and multi-path 
coefficients is non-problematic. The true attenua- 
tion model to be used in this paper is a(p) = p~ 12b 
(see pi). 



5.1 Example 1 

For the purpose of this example, we consider the sen- 
sors to be placed in the first quarter of the x-y plane 
as a spiral array of M = 40 microphones. The spiral 
is represented in a parametric form as 



x M m 




VM m 





4+ ■ 
4 + 



coss 
: sin s 



(33) 



a(p) 



P 1 + E 

e=i 



(32) 



where the angles s are equally spaced in [27r,47r]. 
Our purpose of choosing such sensor arrangement 
was to provide a non-symmetric and still repro- 
ducible arrangement. The sensor locations are 
shown in Fig. [2j The sources used in this exam- 
ple are wideband sources with center frequency 500 
Hz and 200 Hz bandwidth. The sampling frequency 
is 4 KHz. The number of samples available from 
the sources at every sensor is n t = 4000 (i.e., the 
signal duration is 1 second) and the number of fre- 
quency bins is taken to be n/ = 4100. The signal to 
noise ratio (SNR) at every sensor is 20 dB and the 
speed of propagation is considered to be the speed of 
sound as v = 345 m/s. In the proposed minimization 
scheme and more specifically the DE part, we take 
G max = 5. Moreover, we set F = 0.8, Cr = 1 and 
Np = 40. This parameter setting was selected as a 
general DE setting, however more discussions on de- 
termining the DE parameters are available in [19] , 
The general attenuation model is considered to be 
a(p) = jfT 1 + f3ip~ 2 + (32P~ 3 , for which the values (3\ 
and $2 are in charge of tuning the unknown model. 
There is no reverberation in the environment (i.e., 
K = K) and all sensors are synchronized in receiv- 
ing the signal. 

To provide a better understanding of the prob- 
lem, in Fig. [2]the cost function behavior for a known 
attenuation model is shown. In Fig. [2ja the cost 
is shown when the source is located at point (4, 3) 
within the sensors convex hull. All positions are in 
meters. Fig. [2]b shows the cost when the source 
is located at (12,10) outside the sensors region. In 
both cases the cost functions are rather well behaved 
functions away from the sensors. Intuitively, for two 
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Generate an Initial Population 
G = 

Create and Test a New DE 
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Figure 1: A sequential hybridization combing the DE with the LMA 




(a) (b) 

Figure 2: (a) The cost function corresponding to a source located at (4,3) assuming a known attenuation 
model, (b) The cost function corresponding to a source located at (12, 10). 



S 



neighboring points in the domain, sudden variation 
of the cost function with respect to both time delay 
criteria and attenuation model constraints is unlikely 
and hence the resulting cost functions are usually ex- 
pected to be rather slow varying and well behaved 
away from the sensors. 

In Fig. [3] we have shown the iterative procedure 
of finding a single source, once located at (4,3) and 
once at (12, 10). For the first case the source location 
is estimated to be at (4.002,3.008) and the attenua- 
tion coefficients are estimated to be fi\ = -23.85 and 
/?2 = 27.93. In the second case the source estimation 
is (11.999,10.002) and the attenuation coefficients 
are found to be j3\ = 4.19 and /?2 = 1.79. We observe 
that both localization results accurately match the 
exact source positions. The attenuation coefficients 
obtained in both cases are only in charge of fitting 
the low order model to the true model for the source- 
sensor ranges in each problem and due to different 
problems they do not necessarily need to be in the 
same ranges. By providing this low order attenua- 
tion model we provide the flexibility to the problem 
for accurately estimating the sources. 

As a more challenging problem, we consider con- 
current localization of the two sources located at 
(4, 3) and (12, 10). Fig. |4]shows the iterative proce- 
dure of finding the sources. The estimated source lo- 
cations are (4.000, 3.001) and (11.989, 9.993) and the 
attenuation parameters are estimated to be Pi = 4.77 
and fa = -3.94. Again an accurate match between 
the exact source locations and the estimated ones is 
observable. 

We further examine our proposed method 
through a comparison with the AML method devel- 
oped in [13]. For this purpose we start reducing the 
signal samples by reducing the signal duration from 
1 to 0.1 seconds and observing the error caused in 
the source estimation. Here we consider the single 
source localization for the source being located at 
(12,10). Fig. [5] shows the resulting error as the 
signal duration decreases in both methods. As it is 
clearly observed, using both time delay and atten- 
uation information helps our method provide bet- 
ter estimates of the source locations even with less 
available data compared to the AML method which 
only uses the time delay information. We further 
examine the performance of both methods for vari- 
ous SNR values. In Fig. [6] the CRLB is calculated 
for the same single source scenario with the source 
located at (12,10). The RMS errors in estimating 
the x and y components of the source are obtained 



through 50 independent noise realizations for every 
SNR value shown in the figure. Again the proposed 
method shows an acceptable performance regarding 
the closeness to the CRLB and the superior perfor- 
mance compared to the AML method. 



5.2 Example 2 

In a more realistic scenario, we examine the perfor- 
mance of the proposed method in a noisy environ- 
ment where sensor synchronization error and rever- 
beration are likely to happen. The sensor network 
configuration is shown in Figure [7J where three cir- 
cular arrays each composed of 25 sensors centered at 
points (15,5), (2,15) and (5,28) are considered. The 
acoustic source is located at (35,25) and the signal 
specifications are the same as the previous example. 
For this example G max is taken to be 20 to benefit 
more from a global search of a cost function which 
may not be as well-behaved as the previous example 
due to bringing more unknown parameters into the 
problem. The low order attenuation model consid- 
ered in this example is a(p) = p~ l + f3p~ 2 with f3 as 
the tuning parameter. Again an SNR of 20 dB is 
considered at all sensors for all the experiments. 

We first examine the case that the sensors are 
not exactly synchronized to receive the data. For 
this purpose we rewrite the main component of the 
signal in as 

N 



Q^(p m _ n )s n (t T n 

71 = 1 



+ 0, (34) 
where £ is a random variable uniformly distributed 



around zero. Equation (34) basically models the 
asynchronous measurements of the sensor data. In 
Table 1 we have provided the localization results for 
three different synchronization error variances 0.5, 
1 and 2 milliseconds. Clearly the phase error is a 
destructive phenomenon in TDOA localization algo- 
rithms, however, considering the localization errors 
in Table 1, one would observe that exploiting the at- 
tenuation information beside the phase information 
enables our algorithm to perform a rather accurate 
localization task in case of sensors being out of syn- 
chronism. 

Furthermore, a more challenging problem is 
when the reverberation is also taken into account. In 
theory, for the emitted signal to arrive at every mea- 
suring sensor, an individual multi-path filter should 
be considered. Although the formulation in this pa- 
per is general, for the purpose of this example we 
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Figure 3: (a) The progressive estimates of a single source located at (4, 3). The first jumps and good initials 
correspond to applying the DE (b) The progressive estimates of a single source located at (12, 10). (c) The 
evolution of the attenuation model parameters for the single source located at (12, 10) 




(a) (b) 

Figure 4: (a) Concurrent estimates of two sources located at (4,3) and (12,10) (b) The corresponding 
evolution of the attenuation model parameters 
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Figure 5: The localization error verses the signal duration for our proposed method and the AML method 
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Figure 6: (a) The RMS error in estimating the x component of a source located at (12,10). (b) The 
corresponding RMS error in estimating the y component of the source. 
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Figure 7: A sensor network configuration: each array consists of 25 sensors 



have made a reasonable and practical assumption 
that for all the sensors within each array, the filter 
representing the multi-path is identical. In general 
the sensor network may be represented as a collec- 
tion several clusters each composed of sensors closely 
placed and each cluster treated as a single receiving 
node. This assumption prevents dealing with a large 
collection of unknowns (^ m , n ,p an d T m n p ) for every 
source-sensor pair and aggregates them into fewer 
parameters each assigned to the clusters. 

To generate a reverberated signal we use the 
multi-path FIR filters shown in Figure [8] where three 
or four shifted scales of the signal are added to it. 
For the localization purpose, however, we only con- 
sider finding the main indirect path. In other words, 
for every array shown in Figure [7] only one multi- 
path coefficient 7 and one multi-path delay f is to 
be estimated which totally brings 6 unknowns as- 
sociated with the multi-path phenomenon into the 
minimization problem. The remaining minimization 
unknowns are the source coordinates and the atten- 
uation coefficients as before. The fourth row of Ta- 
ble 1 shows the localization result for this problem. 
Although the number of unknowns were relatively 
higher than the previous examples and the cost func- 
tion is clearly not as well-behaved as before, using 
DE as the initial minimization scheme provides a 
suitable starting state for the LMA and this sequen- 
tial technique helps the algorithm make a rather ac- 



curate localization in a noisy and reverberated en- 
vironment. The fifth row of Table 1 corresponds 
to the case of having both the multi-path and the 
synchronization issues, for which the results are still 
promising. The progressive estimates of the target 
throughout the minimization are shown in Figure [9] 



6 Conclusion 

In this paper, we proposed an efficient method for 
localization of multiple wideband sources based on 
both signal attenuation and time delay information. 
The method developed in this paper models the lo- 
calization problem as a minimization problem and 
provides an additional flexibility of not being exactly 
aware of the signal attenuation model. We propose 
a certain function space for the unknown model, and 
tune it iteratively along to estimate the signal source 
locations. The minimization scheme used here is a 
hybrid algorithm, combining the differential evolu- 
tion with the Levenberg-Marquardt algorithm. This 
combination increases the chances of finding a global 
minima while benefits from the speed and computa- 
tional advantages of Newton methods. The accuracy 
and performance of the method is examined through 
several simulations depicting a noisy environment, a 
multi-path environment and lack of synchronization 
among sensors. In the simulations, we compared our 
approach with the approximate maximum likelihood 
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Figure 8: The multi-path impulse response associated with each array. Every multi-path model is assumed 
to hold for all sensors within the corresponding array. The amplitudes are normalized to the amplitude of 
the main signal component. 




(a) (b) 

Figure 9: The progressive estimates of a single source located at (35, 25). The first jumps and good initials 
correspond to applying the DE (a) Considering only the multi-path phenomenon in the environment (b) 
Considering both, the multi-path and sensor synchronization error. 



13 



Table 1: Localization Results for the Sensor Network Configuration in Figure [7] 



Type of Problem 


Synchronization 

Error 
Variance (mS) 


G max 


Number of 

LMA 
Iterations 


Estimated 
Target 
Coordinates 


Localization 
Error 
(meters) 


Array 
Synchronization 


Reverberation 


■ 


□ 


0.5 


20 


21 


(34.62 , 24.91) 


0.386 


■ 


□ 


1.0 


20 


28 


(33.99 , 24.71) 


1.053 


■ 


□ 


2.0 


20 


31 


(33.79 , 24.66) 


1.254 


□ 


■ 





20 


24 


(34.89 , 24.97) 


0.113 


■ 


■ 


0.5 


20 


25 


(34.60 , 24.89) 


0.418 



method which show the superiority of the proposed 
method. 



Using (|37|) in (361 regarding the term 



d(M H M) jdO, would result in 



7 Appendix 

As mentioned earlier, in order to find columns of 
the Jacobian, we are required to find dQ/89, where 
9 is one of the unknown parameters xs n , ys n , Pe, 



dMM^ 



(I-M* H M H ) 



dM 

~de~ 



+ M t fl ^-(/-MMt). 
89 



(38) 



Since Q is a vector containing Based on (35), and knowing (38), we now have 



^Jm^rijp Or T mr) 

sub-vectors Q(f) for / = 0, 1, rif/2, we will only 
find dQ(f)/dO and clearly forming dQ/d9 would be 
aligning the corresponding sub- vectors. 

We first start with replacing the pseudo-inverse 
of K(.f) in (15) which states that 



Q(f) = X(f) 

- K(f)(K h r (/)K(/))- X K fl r (/)X(/). (35) 

We can clearly see that finding dQ(f)/89 re- 

quires having dK(f)(K H (f)k(f))' 1 k" \f)/dO. 
We therefore preliminarily derive some related equa- 
tions. Consider a matrix M, not in general rectan- 
gular, elements of which are dependent on a real 
variable 9. We assume (M H M) exists or in other 

words M f = (M H M)~ M h . Using product rule 
we have 



dMM ] 



d_ 

dO 
dM 



^M(M H M) 1 M H 

. „ ,_j „ d(M H MY 1 „ 
(MM) M h + M^ — 1 M H 



M(M H M) 



-xdM 1 



do 



dM 



M +M 



d{M H M\ 



do 



-M +M 



HdM 1 



(36) 



Also we know that for an invertible matrix M again 
dependent on 9 we have 



dM 

89 



--i3M --i 



(37) 



^tl = (p(f) + P*(f))X(f), (39) 



where 



P(f) = (k iH (f)K H (f)-l)^lk\f). (40) 

To complete the derivation we only need to have 

dK (f)/dQ, which is already discussed in Section 
1X21 
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